#################################################################################
# Replication file for: "A Legislature or a Legislator Like Me? Citizen Demand  #
# for Collective and Dyadic Political Representation"                           #
#                                                                               #
# Jeffrey J. Harden                                                             #
# University of Colorado Boulder                                                #
# jeffrey.harden@colorado.edu                                                   #
#                                                                               #
# Christopher J. Clark                                                          #
# University of North Carolina at Chapel Hill                                   #
# chriclar@email.unc.edu                                                        #
#                                                                               #
# Last update: April 17, 2015                                                   #
#################################################################################
### Packages and data ###
library(foreign)
library(ggplot2)
library(nnet)

## Qualtrics ##
qualtrics <- read.dta("collective-dyadic-qualtrics.dta")

# Treatment conditions #
qualtrics$cydy.r <- ifelse(qualtrics$esblacktreat == "esblacka", 1, 0)
qualtrics$cydn.r <- ifelse(qualtrics$esblacktreat == "esblackb", 1, 0)
qualtrics$cndy.r <- ifelse(qualtrics$esblacktreat == "esblackc", 1, 0)
qualtrics$cndn.r <- ifelse(qualtrics$esblacktreat == "esblackd", 1, 0)

# Black respondents #
q.d <- subset(qualtrics, race == 2)

## CCES ##
cces <- read.dta("collective-dyadic-cces.dta")

# Treatment conditions #
cces$cydy.p <- ifelse(cces$rand_CUB3JH1 == "A", 1, 0)
cces$cydn.p <- ifelse(cces$rand_CUB3JH1 == "B", 1, 0)
cces$cndy.p <- ifelse(cces$rand_CUB3JH1 == "C", 1, 0)
cces$cndn.p <- ifelse(cces$rand_CUB3JH1 == "D", 1, 0)

# Partisans #
cces.d <- subset(cces, as.character(pid3) == "Democrat " | as.character(pid3) == "Republican ")

### Models ###
m.race <- lm(esblacktherm ~ -1 + cydy.r + cydn.r + cndy.r + cndn.r, data = q.d)
m.party <- lm(partytherm ~ -1 + cydy.p + cydn.p + cndy.p + cndn.p, data = cces.d) 
m.party2 <- lm(partytherm ~ -1 + cydy.p + cydn.p + cndy.p + cndn.p, data = subset(cces.d, as.character(race) == "Black"))

### Graphs ###
d.race <- data.frame(trait = "Race (Qualtrics)", pe = coef(m.race), lo = confint(m.race)[ , 1], hi = confint(m.race)[ , 2], cond = c("(A) Collective/Dyadic", "(B) Collective/No Dyadic", "(C) No Collective/Dyadic", "(D) No Collective/No Dyadic"))
d.party <- data.frame(trait = "Party (CCES)", pe = coef(m.party), lo = confint(m.party)[ , 1], hi = confint(m.party)[ , 2], cond = c("(A) Collective/Dyadic", "(B) Collective/No Dyadic", "(C) No Collective/Dyadic", "(D) No Collective/No Dyadic"))
d.party2 <- data.frame(trait = "Party (CCES Blacks)", pe = coef(m.party2), lo = confint(m.party2)[ , 1], hi = confint(m.party2)[ , 2], cond = c("(A) Collective/Dyadic", "(B) Collective/No Dyadic", "(C) No Collective/Dyadic", "(D) No Collective/No Dyadic"))

d.full <- rbind(d.race, d.party)
d.full2 <- rbind(d.race, d.party2)

pdf("levels.pdf")

ggplot(d.full, aes(x = trait, y = pe, fill = cond)) + 
geom_bar(stat = "identity", position = position_dodge(width = .9)) + 
scale_fill_grey(start = .3, end = .7) +
scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, 10)) +
theme(legend.title = element_blank(), legend.position = c(1, 1), legend.justification = c(1, 1), legend.direction = "vertical", legend.text = element_text(size = 15), axis.text = element_text(size = 15), axis.title.y = element_text(size = 18, vjust = 1.5), axis.title.x = element_text(size = 18, vjust = -.1)) +             
geom_errorbar(aes(ymin = lo, ymax = hi), position = position_dodge(width = .9), size = 1, width = .15) +
ylab("Mean Representation Scale Rating") + xlab("Experiment")

dev.off()

pdf("levels2.pdf")

ggplot(d.full2, aes(x = trait, y = pe, fill = cond)) + 
geom_bar(stat = "identity", position = position_dodge(width = .9)) + 
scale_fill_grey(start = .3, end = .7) +
scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, 10)) +
theme(legend.title = element_blank(), legend.position = c(1, 1), legend.justification = c(1, 1), legend.direction = "vertical", legend.text = element_text(size = 15), axis.text = element_text(size = 15), axis.title.y = element_text(size = 18, vjust = 1.5), axis.title.x = element_text(size = 18, vjust = -.1)) +             
geom_errorbar(aes(ymin = lo, ymax = hi), position = position_dodge(width = .9), size = 1, width = .15) +
ylab("Mean Representation Scale Rating") + xlab("Experiment")

dev.off()

### Balance Checks ###
## Qualtrics Data ##
# MNL model #
q.bc <- multinom(as.factor(esblacktreat) ~ factor(votechoice2008) + factor(votechoice2012) + attention2 + party + ideology + econrole + extefficacy + intefficacy + factor(female) + factor(marital) + children + religimp + education + factor(income) + factor(employ), data = q.d)

# Compute p-values #
q.z <- summary(q.bc)$coefficients/summary(q.bc)$standard.errors
q.p <- 2*pnorm(abs(q.z), 0, 1, lower.tail = FALSE)

# Which variable/treatment combination produce a p-value < 0.05? #
q.05 <- which(q.p < .05, arr.ind = TRUE)
cbind(rownames(q.05), colnames(q.p)[q.05[ , 2]])

## CCES Data ##
cces$votechoice2012 <- ifelse(cces$CC354c %in% c("I will not vote in this election", "I'm not sure", "Skipped", "Not Asked"), NA, cces$CC354c) 
cces$faminc2 <- ifelse(cces$faminc %in% c("$500,000 or more", "$150,000 or more", "$250,000 or more", "Skipped"), NA, cces$faminc)
cces$race2 <- ifelse(cces$race %in% c("Native American", "Middle Eastern"), NA, cces$race)
cces$ideo52 <- ifelse(cces$ideo5 == "Not sure", NA, cces$ideo5)
cces$employ2 <- ifelse(cces$employ %in% c("Temporarily laid off", "Unemployed"), NA, cces$employ)
cces$marstat2 <- ifelse(cces$marstat %in% c("Separated", "Domestic partnership"), NA, cces$marstat)
cces$newsint2 <- ifelse(cces$newsint == "Don't know", NA, cces$newsint)
cces$educ2 <- as.numeric(cces$educ)

cces.d <- subset(cces, as.character(pid3) == "Democrat " | as.character(pid3) == "Republican ")

# MNL model #
cces.bc <- multinom(as.factor(rand_CUB3JH1) ~ factor(votechoice2012) + factor(gender) + educ2 + factor(race2) + factor(ideo52) + factor(employ2) + factor(marstat) + factor(pid3) + pew_religimp + factor(child18) + factor(newsint2) + factor(faminc2), data = cces.d)

# Compute p-values #
cces.z <- summary(cces.bc)$coefficients/summary(cces.bc)$standard.errors
cces.p <- 2*pnorm(abs(cces.z), 0, 1, lower.tail = FALSE)

# Which variable/treatment combination produce a p-value < 0.05? #
cces.05 <- which(cces.p < .05, arr.ind = TRUE)
cbind(rownames(cces.05), colnames(cces.p)[cces.05[ , 2]])

### Heterogeneous Treatment Effects? ###
## Race ##
m.race.inc <- lm(esblacktherm ~ factor(esblacktreat) + income + factor(esblacktreat)*income, data = subset(q.d, income < 15))
m.race.att <- lm(esblacktherm ~ factor(esblacktreat) + attention2 + factor(esblacktreat)*attention2, data = q.d)
m.race.edu <- lm(esblacktherm ~ factor(esblacktreat) + education + factor(esblacktreat)*education, data = q.d)
m.race.ecr <- lm(esblacktherm ~ factor(esblacktreat) + econrole + factor(esblacktreat)*econrole, data = q.d)
m.race.ief <- lm(esblacktherm ~ factor(esblacktreat) + intefficacy + factor(esblacktreat)*intefficacy, data = q.d)
m.race.eef <- lm(esblacktherm ~ factor(esblacktreat) + extefficacy + factor(esblacktreat)*extefficacy, data = q.d)

## Party ##
cces.d$faminc3 <- as.numeric(cces.d$faminc)
cces.d$newsint3 <- 5 - abs(as.numeric(cces.d$newsint))
cces.d$educ3 <- as.numeric(cces.d$educ)

m.party.inc <- lm(partytherm ~ factor(rand_CUB3JH1) + faminc3 + factor(rand_CUB3JH1)*faminc3, data = subset(cces.d, faminc3 < 16))
m.party.att <- lm(partytherm ~ factor(rand_CUB3JH1) + newsint3 + factor(rand_CUB3JH1)*newsint3, data = subset(cces.d, newsint3 > 0))
m.party.edu <- lm(partytherm ~ factor(rand_CUB3JH1) + educ3 + factor(rand_CUB3JH1)*educ3, data = cces.d)

# Graphs #
# Attention
x.att <- data.frame(int = 1, rand_CUB3JH1 = rep(c("A", "B", "C", "D"), each = 2), newsint3 = rep(c(2, 4), times = 4))
pred.att <- predict(m.party.att, newdata = x.att, se.fit = TRUE)

d.att <- data.frame(att = rep(c("Only now and then", "Most of the time"), times = 4), pe = pred.att$fit, lo = pred.att$fit - 1.96*pred.att$se.fit, hi = pred.att$fit + 1.96*pred.att$se.fit, cond = rep(c("(A) Collective/Dyadic", "(B) Collective/No Dyadic", "(C) No Collective/Dyadic", "(D) No Collective/No Dyadic"), each = 2))

png("party-att.png", width = 825)

ggplot(d.att, aes(x = cond, y = pe, fill = att)) + 
geom_bar(stat = "identity", position = position_dodge(width = .9)) + 
scale_fill_grey(start = .3, end = .7) +
scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, 10)) +
theme(legend.title = element_blank(), legend.position = c(1, 1), legend.justification = c(1, 1), legend.direction = "vertical", legend.text = element_text(size = 15), axis.text = element_text(size = 15), axis.title.y = element_text(size = 18, vjust = 1.5), axis.title.x = element_text(size = 18, vjust = -.1)) +             
geom_errorbar(aes(ymin = lo, ymax = hi), position = position_dodge(width = .9), size = 1, width = .15) +
ylab("Mean Representation Scale Rating") + xlab("Condition")

dev.off()

set.seed(5234) # Check p-value in treatment D
m <- 1e5
td.lo <- rnorm(m, pred.att$fit[7], pred.att$se.fit[7])
td.hi <- rnorm(m, pred.att$fit[8], pred.att$se.fit[8])
sum(td.lo - td.hi < 0)/m

# Income
x.inc <- data.frame(int = 1, rand_CUB3JH1 = rep(c("A", "B", "C", "D"), each = 2), faminc3 = rep(c(5, 11), times = 4))
pred.inc <- predict(m.party.inc, newdata = x.inc, se.fit = TRUE)

d.inc <- data.frame(inc = rep(c("$40k-$50k", "$120k-$150k"), times = 4), pe = pred.inc$fit, lo = pred.inc$fit - 1.96*pred.inc$se.fit, hi = pred.inc$fit + 1.96*pred.inc$se.fit, cond = rep(c("(A) Collective/Dyadic", "(B) Collective/No Dyadic", "(C) No Collective/Dyadic", "(D) No Collective/No Dyadic"), each = 2))

png("party-inc.png", width = 825)

ggplot(d.inc, aes(x = cond, y = pe, fill = inc)) + 
geom_bar(stat = "identity", position = position_dodge(width = .9)) + 
scale_fill_grey(start = .3, end = .7) +
scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, 10)) +
theme(legend.title = element_blank(), legend.position = c(1, 1), legend.justification = c(1, 1), legend.direction = "vertical", legend.text = element_text(size = 15), axis.text = element_text(size = 15), axis.title.y = element_text(size = 18, vjust = 1.5), axis.title.x = element_text(size = 18, vjust = -.1)) +             
geom_errorbar(aes(ymin = lo, ymax = hi), position = position_dodge(width = .9), size = 1, width = .15) +
ylab("Mean Representation Scale Rating") + xlab("Condition")

dev.off()

set.seed(7957) # Check p-value in treatment D
td.lo <- rnorm(m, pred.inc$fit[7], pred.inc$se.fit[7])
td.hi <- rnorm(m, pred.inc$fit[8], pred.inc$se.fit[8])
sum(td.lo - td.hi < 0)/m

